Cp<-1005 #specific heat capacity
R=287.058 # specific gas constant for dry air (J/(kg·K))
lambda=4
wd<-"C:/Users/marie/OneDrive/Documenten/uhi_max/"
svf<-stack(paste0(wd,"inst/Grids_veg_svf/svf_utrecht_1m.grd"))
fveg<-stack(paste0(wd,"/inst/Grids_veg_svf/greenness_utrecht_smooth_500m.grd"))
fveg<-projectRaster(fveg,crs=crs(svf))
IUtrecht23<-data.frame(lon=52.079,lat=5.139)
coordinates(IUtrecht23)<- ~lat+lon
crs(IUtrecht23)<-crs("+init=epsg:4326")
mapview(fveg) + mapview(IUtrecht23)
## Warning in rasterCheckSize(x, maxpixels = maxpixels): maximum number of pixels for Raster* viewing is 5e+05 ;
## the supplied Raster* has 1348088
## ... decreasing Raster* resolution to 5e+05 pixels
## to view full resolution set 'maxpixels = 1348088 '
IUtrecht23_fveg<-raster::extract(fveg,IUtrecht23)
## Warning in .local(x, y, ...): Transforming SpatialPoints to the CRS of the
## Raster
IUtrecht23_svf<-raster::extract(svf,IUtrecht23)
## Warning in .local(x, y, ...): Transforming SpatialPoints to the CRS of the
## Raster
From the autmaomatic weather stations (AWS) in the rural area the following data is required:
The 10 minute data from global radiation, wind speed, temperature and pressure is used to calculate the UHImax. The hourly relative humidity, precipitation and mean wind are used to filter out synoptic situations with frontal system and fog. Under these different conditions there is no clear relation between the AWS and city temperature.
cabauw_10min<-fread(paste0(wd,"inst/Cabauw_meteo/Cabauw_T_S_W_2016_201809.csv"))
cabauw_10min$IT_DATETIME<-as.POSIXct(cabauw_10min$IT_DATETIME,format="%Y%m%d_%H%M00_000000")
Cabauw.S<-subset(cabauw_10min[which(cabauw_10min$DS_CODE=="348_S_a"),],
select = c("IT_DATETIME","DS_CODE","TOS.Q_GLOB_10"))
Cabauw.T<-subset(cabauw_10min[which(cabauw_10min$DS_CODE=="348_T_a"),],
select = c("IT_DATETIME","DS_CODE","TOT.T_DRYB_10"))
Cabauw.W<-subset(cabauw_10min[which(cabauw_10min$DS_CODE=="348_W_a"),],
select = c("IT_DATETIME","DS_CODE","TOW.FF_SENSOR_10"))
S<-calc_S(Cabauw.S$IT_DATETIME,Cabauw.S$TOS.Q_GLOB_10)
DTR<-calc_DTR(Cabauw.T$IT_DATETIME,Cabauw.T$TOT.T_DRYB_10)
U<-calc_U(Cabauw.W$IT_DATETIME,Cabauw.W$TOW.FF_SENSOR_10)
Meteo_params<-data.frame(DTR,"S"=S$S[1:length(DTR$start)],"U"=U$W)
Meteo_params$start<-as.POSIXct(Meteo_params$start)
Meteo_params$stop<-as.POSIXct(Meteo_params$stop)
saveRDS(Meteo_params,paste0(wd,"inst/Rdata/Meteo_params_cabauw.rds"))
Within the city of Utrecht we have carefully selected three stations from the Wunderground network:
These stations have been measuring the city temperature for the past couple of years. The data can be downloaded using download_time_seq, note that the temperature is in Farenheid and not degrees Celcius. The data has been filtered and interpolated to 10minute timestamps. Unrealistic high and low temperatures are set to NA. In case more than 8 identical measurements in a row occur the value is also set to NA.
cabauw_parms<-readRDS(paste0(wd,"inst/Rdata/Meteo_params_cabauw.rds"))
cabauw_P<-fread(paste0(wd,"inst/Cabauw_meteo/Cabauw_P_10min.csv"))
cabauw_P$IT_DATETIME<-as.POSIXct(cabauw_P$IT_DATETIME,format="%Y%m%d_%H%M00_000000")
cabauw_T<-fread(paste0(wd,"inst/Cabauw_meteo/Cabauw_T_10min.csv"))
cabauw_T$IT_DATETIME<-as.POSIXct(cabauw_T$IT_DATETIME,format="%Y%m%d_%H%M00_000000")
cabauw_PT<-merge(cabauw_P,cabauw_T,by="IT_DATETIME")
air_density<-calc_U(cabauw_PT$IT_DATETIME,(cabauw_PT$TOA.P_NAP_MSL_10*100)/(R*(cabauw_PT$TOT.T_DRYB_10+273.15)))
names(air_density)<-c("start","stop","rho")
cabauw_parms<-merge(cabauw_parms,air_density,by=c("start","stop"))
cabauw_parms$S_new<-cabauw_parms$S/(cabauw_parms$rho*Cp)
cabauw_parms$meteo<-(((cabauw_parms$S_new)*cabauw_parms$DTR^3)/cabauw_parms$U)^(1/lambda)
# svf_grn<-readRDS(paste0(wd,"inst/Rdata/wunderground_svf_grn.rds"))
# df.svf_grn<-data.frame(svf_grn)
# df.svf_grn<-df.svf_grn[df.svf_grn$Station.ID %in% c("IUTRECHT196","IUTRECHT376","IUTRECHT299"),]
Setting all parameters to required to calculate the UHI with the formula of Theeuwes (2017).
# for(i in 1:length(df.svf_grn$Station.ID)){
# STN<-df.svf_grn$Station.ID[i]
STN<-"IUTRECHT23"
stn<-readRDS(paste0(wd,"inst/Wunderground/Filtered/IUTRECHT23_filtered.rds"))
wur_cabauw<-merge(x=stn$interpolated_time,y=cabauw_T,by.x="new_time",by.y="IT_DATETIME")
wur_cabauw<-wur_cabauw[complete.cases(wur_cabauw),]
UHI<-calc_UHImax(time = wur_cabauw$new_time,
Tcity = wur_cabauw$T_int,
Tref = wur_cabauw$TOT.T_DRYB_10)
UHI$Tcity<-calc_U(time=wur_cabauw$new_time,wur_cabauw$T_int)$W
UHI$Tref<-calc_U(time=wur_cabauw$new_time,wur_cabauw$TOT.T_DRYB_10)$W
# UHIcalc_stn<-cbind(cabauw_parms,UHIcalc)
UHI<-merge(UHI,cabauw_parms,by=c("start","stop"))
UHI$svf<-as.numeric(IUtrecht23_svf)
UHI$fveg<-as.numeric(IUtrecht23_fveg)
UHI$Cp<-Cp
UHI$stn<-STN
# p<-ggplot(UHI,aes(UHImeasured,meteo))+geom_point()+geom_abline()
# ggsave(p,filename=paste0(wd,"inst/UHImax/fig/",
# STN,".png"))
write.table(UHI,paste0(wd,"inst/UHImax/",
STN,"_UHIparams.txt"),
row.names = FALSE,
col.names = TRUE,
sep=",")
# }
####################Select control days
meteo.df<-fread(paste0(wd,"inst/Cabauw_meteo/rain_wind_rh_hour.csv"))
meteo.df$IT_DATETIME<-as.POSIXct(meteo.df$IT_DATETIME,format="%Y%m%d_%H%M00_000000")
time<-meteo.df$IT_DATETIME
rh<-meteo.df$BGH.U
rain<-meteo.df$BGH.Q_RH
wind<-meteo.df$BGH.FH
days_subset<-uhi_sub(time=time,wind=wind,rain=rain,rh=rh)
######################
uhimax_files<-list.files(paste0(wd,"inst/UHImax/"),
full.names = TRUE,pattern=".txt")
uhimax_files<-lapply(uhimax_files,fread)
uhimax_Utrecht<-do.call("rbind",uhimax_files)
uhimax_Utrecht<-merge(days_subset,uhimax_Utrecht,by=c("start","stop"))
# uhimax_Utrecht<-uhimax_Utrecht[which(uhimax_Utrecht$Rain==TRUE),]
# uhimax_Utrecht<-uhimax_Utrecht[which(uhimax_Utrecht$Wind==TRUE),]
# uhimax_Utrecht<-uhimax_Utrecht[which(uhimax_Utrecht$rh==TRUE),]
uhimax_Utrecht<-uhimax_Utrecht[which(uhimax_Utrecht$Select==TRUE),]
uhimax_Utrecht<-uhimax_Utrecht[which(uhimax_Utrecht$Tref>17),]
uhimax_Utrecht<-uhimax_Utrecht[which(uhimax_Utrecht$stn %in% c("IUTRECHT23")),]
# uhimax_Utrecht<-uhimax_Utrecht[which(uhimax_Utrecht$stn %in% c("IUTRECHT23","IUTRECHT196","IUTRECHT376","IUTRECHT299")),]
uhi_u<-ggplot(uhimax_Utrecht,aes(U,UHImeasured,colour=DTR)) +geom_point() + scale_color_gradientn(colours = terrain.colors(20))
uhi_DTR<-ggplot(uhimax_Utrecht,aes(DTR,UHImeasured,colour=S_new)) + geom_point()+scale_color_gradientn(colours = heat.colors(20))
uhi_S<-ggplot(uhimax_Utrecht,aes(S_new,UHImeasured,colour=U)) +geom_point() + geom_point()+scale_color_gradientn(colours = topo.colors(20))
ggsave(uhi_u,filename = paste0(wd,"inst/UHImax/fig/uhi_u.png"))
## Saving 7 x 5 in image
ggsave(uhi_DTR,filename = paste0(wd,"inst/UHImax/fig/uhi_DTR.png"))
## Saving 7 x 5 in image
ggsave(uhi_S,filename = paste0(wd,"inst/UHImax/fig/uhi_S.png"))
## Saving 7 x 5 in image
uhi_u
uhi_DTR
uhi_S
ggplot(uhimax_Utrecht,aes(UHImeasured,(2-svf-fveg)*meteo,colour=factor(stn))) +geom_point() +geom_abline() +xlim(0,10)+ylim(0,10)
caret::R2((2-uhimax_Utrecht$svf-uhimax_Utrecht$fveg)*uhimax_Utrecht$meteo,uhimax_Utrecht$UHImeasured)
## [1] 0.2123518
caret::RMSE((2-uhimax_Utrecht$svf-uhimax_Utrecht$fveg)*uhimax_Utrecht$meteo,uhimax_Utrecht$UHImeasured)
## [1] 0.9565242
The formula of Theeuwes(2017) was derived for:
Outside of this range the formula has not been tested.
svf.r<-resample(svf,fveg,method="bilinear")
values(svf.r)[values(svf.r)<0.2 | values(svf.r)>0.95] = NA
values(fveg)[values(fveg)<0 | values(fveg)>0.4] = NA
city_center<-extent(c(5.1,5.13,52.08,52.098))
city_center<-raster(city_center)
values(city_center)<-1
crs(city_center)<-CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
city_center<-projectRaster(city_center,crs=crs(fveg))
st<-(2-svf.r-fveg)*mean(uhimax_Utrecht$meteo)
st<-crop(st,city_center)
mapview(st)